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In conventional molecular simulation, metastable structures often survive over considerable computational time, result- 
ing in difficulties in simulating equilibrium states. In order to overcome this difficulty, here we propose a newly devised 
method, molecular Monte Carlo simulation of systems connected to three reservoirs: chemical potential, pressure, 
and temperature. Gibbs-Duhem equation thermodynamically limits the number of reservoirs to 2 for single compo- 
nent systems. However, in conventional simulations utilizing 2 or fewer reservoirs, the system tends to be trapped in 
metastable states. Even if the system is allowed to escape from such metastable states in conventional simulations, the 
fixed system size and/or the fixed number of particles result in creation of defects in ordered structures. This situation 
breaks global anisotropy of ordered structures and forces the periodicity of the structure to be commensurate to the 
system size. Here we connect the such three reservoirs to overcome these difficulties. A method of adjusting the three 
reservoirs and obtaining thermodynamically stable states is also designed, based on Gibbs-Duhem equation. Unlike 
the other conventional simulation techniques utilizing no more than 2 reservoirs, our method allows the system itself to 
simultaneously tune the system size and the number of particles to periodicity and anisotropy of ordered structures. Our 
method requires fewer efforts for preliminary simulations prior to production runs, compared with the other advanced 
simulation techniques such as multicanonical method. A free energy measurement method, suitable for the system with 
the three reservoirs, is also discussed, based on Euler equation of thermodynamics. This measurement method needs 
fewer computational efforts than other free energy measurement methods do. 



I. INTRODUCTION 

Metastable structures are found in a variety of physical 
systems, e.g. glasses, amorphous solids of colloids 1 ' 2 , stalk 
intermediate structures in biological membrane fusion pro- 
cess 3 ~ 6 , and so on. Metastable states correspond to regions 
of the local free energy minima in phase space. In conven- 
tional simulation of the canonical ensemble, once the sys- 
tem is captured in these regions, the system often remain in 
these non-equilibrium states for enormously long time. Such 
metastable states are frequently found in macromolecular and 
colloidal systems, and make the simulation studies on equi- 
librium states difficult. Even if the system can escape from 
the non-equilibrium states to the ordered equilibrium state, 
the constant number of particles and the constant system size 
cause defects in the ordered structure. This situation breaks, 
in a global scale, the anisotropy of the ordered structure and 
forces the periodicity of the ordered structure to be commen- 
surate to the system size. In order to find defect-free equi- 
librated ordered structures, we should finely tune the system 
box size (l x , L y , L 2 ) as well as the number of particles N so 
that both the anisotropy and the periodicity of the ordered 
structure, which are not known a priori, are not violated by 
the periodic boundary conditions. This fine tuning for the or- 
dered structure is in general a tedious task. 

Advanced simulation techniques, e.g. multicanonical en- 
semble method 7 ' 8 and umbrella sampling 8-10 , allow us to al- 
most homogeneously sample the whole phase space at con- 
stant N and \L X , L y , L 2 ), with the help of artificial weights that 
reduce the occurrence probability of non-equilibrium states. 
However, microscopic states in equilibrium, sampled by these 
advanced techniques, are restricted to microstates with the 
given set of constant N and constant (l x ,L v ,L z ). Free en- 



ergy landscapes of the system with different sets of N and 
(L x ,L v ,L z ^j at the same particle density are not searched. 
These extensive variables simultaneously need manual fine 
tuning for the purpose of finding the most stable state in these 
free energy landscapes. For example, both N and (l x , L y , L z ) 
of perfect crystals should be integer multiples of the unit 
structure. Furthermore these advanced methods also require 
an advanced programming and large amounts of complicated 
preparation, e.g. accurate calculation of free energy 11 and the 
precise adjustments of the artificial weights, prior to the pro- 
duction simulation runs. In addition, such unphysical sam- 
pling processes with the artificial weights make it difficult to 
trace physical trajectories in the phase space. 

Here we devise molecular Monte Carlo simulation method 
of systems connected to three reservoirs (hereafter we call it 
"three-reservoirs method") 12 , chemical potential /i, pressure 
P, and temperature T, for seeking the most stable states of the 
target systems, i.e. the equilibrium structures. Due to Gibbs- 
Duhem equation, 

S dT -VdP + Ndp = 0, (1) 

where S is the entropy and V is the volume, the number of 
reservoirs is thermodynamically limited to no more than 2 for 
single component systems. However, we connect the three 
reservoirs in order to overcome the above difficulties of the 
other conventional and advanced simulation techniques. In 
order to perform this, we introduce a method for adjusting 
these three reservoirs and obtaining thermodynamically stable 
states. 

The total number of particles N and the system box size 
(L x ,Ly,L^} are additional degrees of freedom of the system 
connected to these three reservoirs. These additional degrees 
of freedom correspond to additional dimensions of the phase 
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space, which provide shortcuts from the non-equilibrium state 
to the equilibrium state. In addition, unlike the other simula- 
tion techniques utilizing 2 or fewer reservoirs, these degrees of 
freedom allow the system itself to simultaneously tune N and 
(L x ,L y ,L^, so that the system reaches the true equilibrium 
ordered structure. Furthermore, our method requires fewer 
efforts for the preliminary simulation prior to the production 
simulation runs. 

Guggenheim formally introduced Boltzmann factor (sta- 
tistical weight) of the ensemble with the three reservoirs 13 . 
Prigogine and Hill also studied the same ensemble later 14,15 . 
These early works, however, focused on mathematical aspects 
of the partition function, i.e. mathematical formalism of the 
ensemble, since their goal was to discover a universal and 
generalized expression for a partition function applicable to 
any thermodynamically acceptable ensembles 16 ' 17 . By con- 
trast, physical aspects of the ensemble were wholly left for 
the future. In the present work, we study the physical as- 
pects of this ensemble intuitively and thought-experimentally. 
In addition, we also analytically solve maximization prob- 
lems of entropy densities, which corroborates our intuitive 
and thought-experimental study. Finally, we design three- 
reservoirs method based on these physical aspects of the en- 
semble and show simulation results on non-trivial globally- 
anisotropic defect-free ordered structures of colloidal sys- 
tems. 

In spite of Gibbs-Duhem equation, a system connected to 
the three reservoirs can be realized in experiments. For exam- 
ple, we can imagine a system that obeys the grand canonical 
ensemble (i.e. constant //, V, and T), and replace one of its 
walls with a free piston facing to a reservoir of pressure P, i.e. 
the 3rd reservoir. In the present article, we theoretically con- 
struct thermodynamics and statistics of the system connected 
to the three reservoirs. 

Based on Euler equation of thermodynamics, we also pro- 
pose a method for measuring entropy and free energy directly 
from the simulations of the systems connected to the three 
reservoirs. This measurement needs fewer computational ef- 
forts than the other free energy calculation methods using 
molecular simulation. 

We design the algorithms of the three-reservoirs method 
based on conventional Monte Carlo (MC) simulation meth- 
ods of the grand canonical ensemble (juVr-ensemble) and the 
isothermal-isobaric ensemble (A^Pr-ensemble). We give a 
brief description of these conventional molecular MC tech- 
niques in appendix A. Thermodynamics, statistical mechan- 
ics, and simulation methods of the system with the three reser- 
voirs are studied in section II. Finally, we summarize the 
present work in section III. 



II. //Pr-ENSEMBLE 

Here we discuss the basic formalism of the three-reservoirs 
method, where the method for adjusting the three reservoirs is 
also developed. Thermodynamic properties of this system are 
discussed in section II A. A microscopic formulation of three- 
reservoirs method based on statistical mechanics will be given 



in section II B, where we solve the maximization problem of 
the statistical entropy per volume. Algorithms for the simu- 
lation based on this statistical formulation are constructed in 
section II D. Simulation results to demonstrate efficiency and 
stability of three-reservoirs method are given in section HE. 
Finally, a new and simple technique to measure the entropy 
and the free energy of the system is proposed in section II G. 



A. Thermodynamics in yu/T-ensemble 

As an example of thermodynamic systems, we consider a 
gas contained in a diathermal box with a free piston. This 
box is placed in an environment at constant T and constant 
P. These two intensive variables (T, P) determine the other 
intensive variables of this system, e.g. the chemical poten- 
tial //, the number density of particles p = N/V, and the free 
energy per particle G/N. This means that thermodynamic de- 
grees of freedom of this system are equal to 2, which results 
from Gibbs-Duhem equation, eq. (1). In conventional simula- 
tion methods, N is also fixed at some value (A^Pr-ensemble), 
whereas states and phases of the system are independent of N, 
i.e. a change of N only scales the extensive variables of the 
system. In other words, phase diagrams constructed in PT- 
plane are independent of the extensive variables. Instead of 
fixing this insignificant N, we connect this system to a reser- 
voir of //, whose value is determined by (T, P) through Gibbs- 
Duhem equation, eq. (1). Since Gibbs-Duhem equation is sat- 
isfied, this third reservoir does not affect the equilibrium state 
of the system. This condition corresponds to a thermodynam- 
ically stable point based on Gibbs-Duhem equation, which 
results in an equation of state that links T, P, and /i. At a 
thermodynamically stable point of this ensemble, the exten- 
sive variables of the system, e.g. V and N, are freely scaled, 
i.e. indeterminate and fluctuating, while the system keeps all 
the intensive variables fixed. In simulating systems connected 
to the three reservoirs at the thermodynamically stable points, 
we can choose simulation runs at small N, which are computa- 
tionally advantageous. Gibbs free energy per particle, which 
must be minimized in A^Pr-ensemble before adding the 3rd 
reservoir of constant ,u, is unchanged even after this 3rd reser- 
voir is connected to the system. 

In the above example, (T, P) is given from the outside of 
the system; /j. is adjusted according to (T, P) and connected 
to the system as an additional reservoir. Two other combi- 
nations (T,jj) and P, and (jj., P) and T also work in a similar 
manner. Therefore, in addition to Gibbs free energy per parti- 
cle, both grand potential per volume and the thermodynamic 
potential of [iPS -ensemble per S , i.e. (E-/iN+PV)/S , are si- 
multaneously minimized in the system connected to the three 
reservoirs, where S and E are entropy and internal energy re- 
spectively. 

Gibbs-Duhem equation explains this simultaneous mini- 
mization of the 3 free energy densities. When we add the 3rd 
reservoir to the system and set the system at the thermody- 
namically stable point, the intensive variable of the 3rd reser- 
voir needs to be adjusted, in advance of the connection of this 
additional reservoir. Through this adjustment of the 3rd inten- 



3 



1 1 




P 


O 


* V,N 


r 


1 *^ 


■° * 





T 



FIG. 1. Sketch of a system connected to the three reservoirs. The 
reservoir 1 and the system are allowed to exchange particles to fix 
the chemical potential of the system at fi, where N of the system is 
a dynamic variable. Between the reservoir 2 and the system, a free 
piston is placed. This piston moves and changes the volume V to fix 
the pressure of the system at P. The system and these reservoirs are 
connected to a thermostat at T, the reservoir 3. 



sive variable, the corresponding free energy density is mini- 
mized. The same system at the same thermodynamically sta- 
ble point is also constructed by the two other combinations of 
the 3 intensive variables. In other words, NPT, fiVT, and 
pPS -ensembles simultaneously underlie the ensemble with 
the three reservoirs. This results in the simultaneous mini- 
mization of the 3 free energy densities. In the other ensembles, 
however, any sets of corresponding 3 external parameters, e.g. 
(T, V,N) in the canonical ensemble (A^Vr-ensemble), can be 
selected arbitrarily and the adjustment of the external param- 
eters is not demanded. Therefore, the corresponding free en- 
ergy, e.g. Helmholtz free energy in A^VT-ensemble, is mini- 
mized in a conventional ensemble, whereas other free energies 
are not minimized due to the absence of underlying ensem- 
bles. 

A system connected to the three reservoirs is sketched in 
Fig. 1. The reservoirs 1 and 2 define the values of p and P 
of the system respectively. The particle density of these 2 
reservoirs, denoted by p^J s and p® respectively, determine the 
p and P. The system and the reservoirs 1 and 2 are connected 
to a thermostat at T, the reservoir 3. In appendix A, we give 
a derivation of this p, and its explicit expression is given in 
eq. (Al). 

As another simple example, we thermodynamically con- 
sider a system composed of a single-component ideal gas. 
In this example, we assume that the reservoirs are also com- 
posed of the same ideal gas. At the thermodynamically stable 
point, a relation p[.'s = pill holds and the particle number den- 
sity of the system, p, also equals these particle number den- 
sities of the reservoirs; i.e. p[H - p® = p. However, when 
Pies > Pres' both N and V diverge, since the reservoir 1 con- 
tinues to increase N of the system aiming for a large p value 
and the reservoir 2 increases V aiming for a small value of p. 
When p'es < p[es, both N and V vanish. Therefore, the sys- 
tem reaches equilibrium only at the thermodynamically stable 



point. In other words, outside the thermodynamically stable 
point, the system is always in non-equilibrium and both the 
intensive and the extensive variables are indeterminate. These 
results also apply to systems composed of interacting particles 
(e.g. non-ideal gases). We utilize these divergence and vanish- 
ment of the target system as a criterion for the equilibration, 
which can be used for the automatic adjustment of the three in- 
tensive variables, e.g. by a bisection method, to determine the 
thermodynamically stable point. The system quickly diverges 
or vanishes outside the vicinity of the thermodynamically sta- 
ble points 12 . The speed of the divergence and the vanishment 
increases with the difference in the intensive parameter sets 
from the thermodynamically stable point. 

On the other hand, when we need a long simulation run in 
the vicinity of the thermodynamically stable point, indetermi- 
nate N could cause a computational problem, since the exten- 
sive variables could become extremely large or vanish. How- 
ever, an appropriate choice of N and V makes the system last 
for a long time, within which good statistics of simulation re- 
sults, e.g. particle density and lattice constants of crystals, are 
obtained 12 . We can determine such simulation results with ac- 
curacy enough to obtain the equilibrated structure of the sys- 
tem with the three reservoirs. If we need a far longer simula- 
tion run, the ensemble can be switched to one of the conven- 
tional methods, e.g. NPT-ensemble or A^VT-ensemble. As 
these switched ensembles are free from the problem of the in- 
determinate N, they allow us to perform a longer simulation 
run of the equilibrated structure obtained via three-reservoirs 
method. 

In the present article, we tentatively call the ensemble of 
the systems connected to the three reservoirs yuPr-ensemble, 
since this is obtained as an equilibrium condition between the 
three intensive variables and as a combination of pVT, NPT, 
and pPS -ensembles. 



B. Statistical mechanical properties of particles in 
//Pr-ensemble 

Here the statistical mechanical properties of particles in 
//Pr-ensemble are discussed. According to Gibbs-Duhem 
equation, the thermodynamic potential of pPT -ensemble is 
identically equivalent to zero in the thermodynamic limit, 
whereas systems obeying this ensemble have certain degrees 
of freedom in the phase space. This means that Boltzmann 
factor fluctuates in statistical mechanics and that the parti- 
tion function of this ensemble, i.e. summation of statistical 
weights over the whole phase space, is defined. This Boltz- 
mann factor is calculated in the present section as a natu- 
ral extension of those for pVT and A^Pr-ensembles, and de- 
termines detailed balance conditions necessary for designing 
three-reservoirs simulation method. 

For determining this Boltzmann factor, we solve a max- 
imization problem of the statistical entropy per volume in 
/ifT-ensemble. The statistical entropy is defined as 16 ' 18 , 

-kn^jPjlogpj, (2) 

j 
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where the suffix, j, represents microstates of the system, 
pj denotes corresponding occurrence probability, and k B is 
Boltzmann constant. Because the extensive variables of pPT- 
ensemble is indeterminate, we choose the entropy density to 
be maximized. With the use of eq. (2), the statistical entropy 
per volume is given by, 



PjlogPj 
Vj ' 



(3) 



The equilibrium probability distribution, {p\,p2, ■ ■ • ,Pj, ■ ■ ■ }, 
which maximizes eq. (3) under the constraints in pPT- 
ensemble, is determined by: 
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(4) 
(5) 



where (• • • ) is the thermal average specified according to the 
reservoirs. The constraint of eq. (4) represents the normaliza- 
tion condition. The two constraints given by eq. (5) instead 
of three constraints are due to the thermodynamic degrees of 
freedom, which equal 2. With the use of Lagrange multipliers, 
a, ft, and y, this maximization problem is reduced to 18 , 
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(6) 



This reduced problem, eq. (6), is solved by a partial derivative 
with respect to pj. The solution is: 



P j = Gxp[aVj-fiEj+yNj-i\, 



(7) 



where a = a/k B , /3 = -fi/k B , and y = y/k B . Equations (3) - 
(5), and (7) give the statistical entropy per volume determined 
as thermal average as: 



S 
V 



(8) 



Equating eq. (8) with the thermodynamic relation for the ther- 
modynamic potential of /L/Pr-ensemble, denoted by <£, 



<S) = E-TS -fiN + PV, 



we obtain: 



a - - 



k B T' 



P = 



1 



k B T 
k B T log e. 
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k B T ' 



(9) 

(10) 
(11) 



In the thermodynamic limit, the right-hand side of eq. (1 1) is 
essentially zero compared with the other extensive variables in 
eq. (9). This result coincides with the Euler equation in ther- 
modynamics. Furthermore, for finite V, the last term of eq. (8) 



is monotonically increasing with decreasing V. This indicates 
that due to the principle of maximizing entropy a finite size 
system as is used in the computer simulation has a tendency 
to shrink even at the thermodynamically stable point. Finally, 
from eq. (7), the probability distributions of the microstates 
are obtained as, 



1 

pj = - exp 



k B T 



-Vi- 



1 

k^f 1 



k B r 



■Ni 



(12) 



We confirmed that these results are also obtained by maximiz- 
ing the statistical entropy per particle or the statistical entropy 
per internal energy instead of the statistical entropy per vol- 
ume as has been done in eq. (6). 

The Boltzmann factor of /iPr-ensemble is, 



exp 



1 

'k^f 



{PV-pN+U(r u ..., r N ) + E kinetic } 



(13) 



where ^kinetic denotes the total kinetic energy of the system, 
U the potential energy of the system, and r, is the spatial co- 
ordinates of the particle ;. As we have seen above, this is a re- 
sult of a natural extension of pVT and A^Pr-ensembles. This 
Boltzmann factor determines the statistical properties of sys- 
tems equilibrated with the three reservoirs and is equivalent 
to the Boltzmann factor of A^Pr-ensemble at fixed N, and is 
equivalent to the Boltzmann factor of //VF-ensemble at fixed 
V. 



C. Thermodynamic potential of yu/T-ensemble 

Guggenheim formally derived the Boltzmann factors (sta- 
tistical weight) of various ensembles 13 . Assuming that the en- 
semble averages of extensive variables, e.g. V and N, were 
determined in yuPr-ensemble, Guggenheim also introduced 
statistical weight of this ensemble, based on analogy between 
other conventional ensembles 1317 . However, in the calcula- 
tion of the maximization of entropy density discussed in sec- 
tion II B, Guggenheim's assumption corresponds to keeping 
the averages (E) and (N) fixed, instead of the constraints 
eq. (5). Guggenheim's assumption contradicts the indetermi- 
nation of the extensive variables, as was pointed out by Pri- 
gogine and further discussed by Sack 1416 . Prigogine showed 
that the summation of the resulting Boltzmann factor over 
the phase space diverges and therefore concluded that the 
resulting partition function of //PF-ensemble does not have 
any physical meanings 1416 . The thermodynamic potential of 
/iPr-ensemble, <D, determined from such partition function 
could be indefinite while it should identically equal zero in the 
thermodynamic limit. The true thermodynamic potential that 
dominates this ensemble, similar to the Helmholtz free energy 
in A^Vr-ensemble, remains unknown since Guggenheim's ar- 
ticle. In the following, we will answer Prigogine's criticism 
and give an explicit expression of the thermodynamic poten- 
tial for //Pr-ensemble. 

In conventional ensembles, the statistical weight takes non- 
zero values only in the vicinity of the averages (N) or (V) in 
the phase space. Outside this vicinity, the statistical weight 
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quickly decreases to zero. This suppresses the divergence of 
the partition functions, i.e. the summation of the statistical 
weights, in the case of conventional ensembles. On the other 
hand, in //PF-ensemble, there is no such limitation because of 
the indeterminate extensive variables. The statistical weight 
of //Pr-ensemble at each microstate keeps non-zero values at 
any N or V. This results in the divergence of the partition 
function of /iPF-ensemble. However, ratios of the statistical 
weights between any pair of microstates are still defined. This 
feature guarantees the physical validity of //PF-ensemble. In 
this case, the trajectory of the system in the phase space is 
similar to a free random walk in infinitely large space without 
boundaries. 

Moreover, in the present study, equations (11) and (12) 
calculated with the constraints eqs. (4) and (5) indicate that 
the summation of the statistical weights equals e for pPT- 
ensemble. Therefore, the corresponding thermodynamic po- 
tential <£, eq. (1 1), is negligibly small compared with the other 
extensive variables in eq. (9) and vanishes in the thermody- 
namic limit. Furthermore, by the thermodynamic considera- 
tion given in section II A, we have shown that /iPT -ensemble 
is obtained by combining 3 underlying ensembles each with 
2 reservoirs, i.e. NPT, pVT, and pPS -ensembles. This ther- 
modynamic consideration means that free energy densities of 
these underlying ensembles, i.e. Gibbs free energy per parti- 
cle, grand potential per volume, and (E - pN + PV)/S, are 
simultaneously minimized in yuPr-ensemble, rather than <£>. 
This corresponds to the minimization of Helmholtz free en- 
ergy in A^Vr-ensemble. 

D. MC simulation method in \iPJ -ensemble 

The present simulation method is constructed based on 
conventional simulation methods in the grand canonical en- 
semble (//Vr-ensemble) and jVPr-ensemble. The simula- 
tion algorithms of the particle insertion and deletion in pVT- 
ensemble (see appendix A 1) and the system size change in 
Wr-ensemble (see appendix A 2) are directly utilized in our 
method. This compatibility between the present method and 
the conventional methods demonstrates that the algorithms of 
the present method satisfy detailed balance condition. 

One simulation step of the present method is composed of 
the following 4 trial steps: 

i) with probability pa/2, trial particle insertion into the sys- 
tem, 

ii) with probability pa/2, trial particle deletion from the sys- 
tem, 

iii) with probability p v , trial system size change, 

iv) with probability 1 - po - pv, trial displacement by 
Metropolis algorithm i.e. perturbation of one particle, 

is chosen, where p c and p v are constants fixed in an interval 
< pcpv < 1. 

With the use of the simulation algorithms of trial particle in- 
sertion and deletion in /iVF-ensemble, the insertion and dele- 
tion of the present method, steps i) and ii), are performed. 



During this particle exchange between the system and the 
reservoir 1, the system size (L X ,L V ,L Z ) is fixed. This particle 
exchange in the present method satisfies the detailed balance 
condition because it is guaranteed in /iVF-ensemble. 

The trial system size change, step iii), is performed with use 
of the simulation algorithms in iVPr-ensemble, during which 
N is fixed. This system size change in the present method also 
satisfies the detailed balance condition. Unlike the conven- 
tional MC simulations in A^Pr-ensemble based on McDon- 
ald's method 10 , L x , L y , and L z are independently changed in 
the present method. 

The trial move of one particle, step iv), is performed by 
Metropolis algorithm in jVVr-ensemble, which also satisfies 
the detailed balance condition. 

Therefore, the present MC simulation method for pPT- 
ensemble fulfills the principle of detailed balance. See also 
section II D 1 . 

Our algorithm indicates that a short-time average of an in- 
tensive physical quantity in //PF-ensemble is approximated 
by ensemble averages of pVT and Wr-ensembles at corre- 
sponding N and V, which is discussed in appendix B. 

Our simulation is performed in a rectangular system box 
with independently changing system size (L x , L y , L z ). Since 
any crystals fit into rectangular boxes with periodic boundary 
conditions, we do not have to introduce Parrinello-Rahman 
method 1019 21 , which allows the change in the shape of the 
simulation box. 

Our simulation algorithm is similar to Gibbs ensemble tech- 
nique, which is utilized for simulation of phase equilibria in 
A^Vr-ensemble 8 ' 10 ' 22 ' 23 , where the phases coexisting in the 
same system box in A^Vr-ensemble exchange both volume 
and particles. A system connected to the three reservoirs at a 
thermodynamically stable point corresponds to this Gibbs en- 
semble, when one of the coexisting phases in Gibbs ensemble 
is assumed to be infinitely large that plays the role of the two 
reservoirs p and P. 



1. Detailed balance condition and ergodicity of 
three-reservoirs method. 

In this section, the detailed balance condition and the er- 
godicity of three-reservoirs method are discussed. Steps i) 
and ii) change N according to the detailed balance condition 
as in the same way that the yuVr-ensemble does. Step iii) 
changes (L x , L y , L z ) according to the detailed balance condi- 
tion as in the A^Pr-ensemble. The particle coordinates are 
updated in step iv), by the standard Metropolis algorithm of 
jVVr-ensemble. As a result, three-reservoirs method satis- 
fies the principle of detailed balance and ergodicity, based on 
the conventional ensembles which fulfill ergodicity. This in- 
dicates that N, (l x , L y , L z ), and the particle coordinates are 
simultaneously updated in the phase space, so that the sys- 
tem realizes the equilibrium state. This also means that statis- 
tics of //Pr-ensemble contradicts none of the three underly- 
ing ensembles with 2 reservoirs, i.e. pVT, NPT, and pPS- 
ensembles. 
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E. Examination of three-reservoirs method 

In this section, we demonstrate the efficiency and stability 
of the three-reservoirs method using several examples. 

The first example is a model polymer-grafted colloidal sys- 
tem, which was observed to show various exotic metastable 
phases each of which has a long life time 2 . 

Colloidal particles are made from metals, polymers, etc. 
and are often modeled with hard spheres 24 . On the other 
hand, owing to van der Waals attraction acting on particle 
surfaces, colloids are aggregating and make precipitates af- 
ter a long time. For the purpose of stabilizing colloidal dis- 
persions against the precipitation, linear polymer chains are 
often grafted onto the surfaces of the colloids. These are 
called polymer-grafted colloids. Depending on the physi- 
cal and chemical properties of the grafted polymers, interac- 
tion between polymer-grafted colloids changes significantly 1 . 
Polymer-grafted colloids have several industrial applications 
due to such useful characteristics, for example filler particles 
immersed in a polymer matrix, and the particles in electro and 
magnet rheological fluids. 

In our previous work 2 , we studied the phase behavior of 
colloidal particles onto which diblock copolymers are grafted. 
Pair interaction potential between these polymer-grafted col- 
loids was numerically determined via self-consistent field cal- 
culation 2 as a function of the distance between centers of the 
particles, r. This potential has been approximated by spheri- 
cally symmetrical repulsive square-step potential with a rigid 
core of diameter <t\ and a square-step repulsive potential of 
diameter cr 2 and height e as: 

(f>(r) - co r < <T\, 

<f>(r) = £ Q (> 0) <T\ < r < <x 2 , 
(p{r) = cr 2 < r. 

where the step potential is originated from the grafted polymer 
brushes. This interaction potential <p(r) is purely repulsive. 
Simulating particles interacting via (p(r) in /VVr-ensemble, 
we have studied phase behavior of these colloidal systems. 
These MC simulation results 2 show that, at low temperature, 
high pressure, and cr 2 /cr l ~ 2, our particles self-assemble into 
string-like assembly. The positions and the mean-square dis- 
placement of the particles show that the string-like assembly 
is observed in disordered solid phases. Actually, such a string- 
like assembly has been observed recently in experiments 25 . In 
addition to such string-like assembly, various structures, e.g. 
dimers and lamellae 26 ' 27 , and also glass transition 28 are ob- 
served in the same model system. It was also shown that the 
particles interacting via continuous repulsive potential simi- 
lar to the above <f>(r) show these string-like and other various 
assemblies 29 . In these recent studies using /VVr-ensemble at 
finite T in both 2 and 3-dimensions, it was shown that this 
string-like assembly with a local alignment in the same di- 
rection but with a global isotropy is a metastable structure. 
Although a variety of ground states of the same model sys- 
tem at zero temperature have been discovered via genetic al- 
gorithms 30 , equilibrium states at finite temperature have not 
been understood yet. 



We simulate the equilibrated states of our model system at 
finite T via three-reservoirs method. In the present simulation 
work, cti and eo are taken as the unit length and the unit en- 
ergy, respectively. Simulation is performed on 2 dimensional 
systems. We define dimensionless chemical potential as 6 , 

p' :=fi/k B T-log(A 2 loi). (14) 

The thermal de Broglie wave length A, defined in eq. (A3) 
in appendix A, is removed from p in eq. (14), since simula- 
tion results are independent of A, which is discussed in ap- 
pendix A. In the present article, o^/cn = 2.0 is fixed. 

In the initial state, No particles are arranged on a homoge- 
neous triangular lattice in a square system box, L x /L y = 1.00, 
with the periodic boundary condition. For the trial move of 
the particles, i.e. Metropolis algorithm, a particle is picked at 
random and given a uniform random and isotropic trial dis- 
placement within a square whose sides have length 0.4o-!. AL 
is fixed at O.Olcr,. The probability p c = 0.1 and p v = l/N () . 
With this pv, the computational time for the simulation is 
about twice as long as the simulation in /VVr-ensemble. 1 
Monte Carlo step (MCS) is defined as N simulation steps. 
The Mersenne Twister algorithm 3133 is adopted as a uniform 
random number generator for our simulation. 



1. String-like assembly. 

In our previous /VVr-ensemble simulation with the poten- 
tial step width, o- 2 /cr 1 = 2.0, we found 2 that the string length 
diverges at ksT/eo = 0.12 and per 2 « 0.451. At this low 
ksT leo, using three-reservoirs method, we simulate the sys- 
tem. In the present simulation, the density per 2 is initially set 
at 0.451 for/Vo = 1254 system. 

First, simulating the system at various values of p.', 
we search for the thermodynamically stable point at fixed 
Pcr 2 /e = 1.0. The given//' is, 1): p' = 27.63, 2): p' = 29.93, 
3): p' = 31.54, 4): p' = 32.23, and 5): p' = 34.53. 
As an example, a snapshot of the system at 1) is presented 
in Fig. 2(a). Despite different p', all the systems of 1) to 
5) show similar well-aligned globally-anisotropic defect-free 
string-like assembly, though the system in /VVr-ensemble at 
pu\ = 0.451 and N = 1200 presents the globally-isotropic 
string-like assembly as is shown in Fig. 2(c) 2 . Only small 
and short-lived defects caused by the thermal fluctuation can 
be generated in the systems simulated with yuPr-ensemble, 
whereas many long-lived defects are observed in NVT sim- 
ulation. Time evolutions of per 2 and TV at 1) to 5) are given 
in Fig. 3. All the data for the time evolution of pcr\ shown in 
Fig. 3(a) are fluctuating in the vicinity of pcr\ ~ 0.451, regard- 
less of p' . It would be worth noting that the simulations started 
from different initial conditions, e.g. different initial particle 
density per 2 and different system aspect ratio of the simula- 
tion box L x /L y , also reach the same results as long as the in- 
tensive variables of the reservoirs, p', Pcr 2 /e , and ksT/eo are 
the same. Furthermore, outside this range of p! from 1) to 5), 
the system diverges (TV — > oo) or vanishes (TV — » 0) just after 
the simulation starts. These results illustrate that this region 
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(a) (b) (c) 

FIG. 2. Snapshots of the system simulated in yiiPr-ensemble. (a), k B T/eo = 0.12, Po^/eo = 1.0 and p' = 27.63. (b), k B T/eo = 0.12, 
Po"i/eo = 1.1 and p! = 29.93. Black dots represent the centers of the particles and grey lines denote networks of overlaps between the 
particles. A snapshot of the system simulated in iVVr-ensemble at k B T/eo =0.12, pa\ = 0.451, and A' = 1200 is also shown in (c) 2 . 
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FIG. 3. Time evolution of p = N/V, (a), and time evolution of ./V, 
(b), simulated at k B T/e = 0.12 and various {Pa\le , p'). Lines 1) 
to 5) are simulation results at fo" 2 /e = 1-0. 1): p' = 27.63. 2): 
p' = 29.93. 3): p' = 31.54. 4): p' = 32.23. 5): p' = 34.53. Lines 
6) to 8) are simulation results at p' = 29.93. 6): Pa-]/e Q = 0.95. 7): 
Pa- 2 Je = 0.98. 8): />o- 2 /e = 1.1. 



of p! is located in the vicinity of the thermodynamically stable 
point at this pressure, Per 2 /eo = 1 .0. 

On the other hand, time evolution of N, plotted in Fig. 3(b), 
indicates that the total system size depends on p! in a system- 
atic manner. At small MCS, line 1) indicates the tendency of 
the vanishment and lines 2) to 5) the tendency of the diver- 
gence. This means that, in a short computational time, the 
thermodynamically stable point is expected to lie between 1): 
p! = 27.63 and 2): p' = 29.93, which corresponds to a relative 
error of some percent. The abrupt time evolution of N stops at 
less than 10 7 MCS, whereas the system slowly continues di- 
verging or vanishing at larger MCS. The MCS needed for the 
divergence or the vanishment becomes larger when we choose 
parameter sets close to the exact values at the thermodynami- 
cally stable point. This is utilized as a criterion for measuring 
the convergence of the thermodynamic intensive variables in 
//Pr-ensemble. Although we stop the simulation with this 
relative error of some percent, the accuracy of the thermody- 
namically stable point could be raised, e.g. by the bisection 
method improving the accuracy of p' . 

Lines 6) to 8) in Fig. 3 show the results of a similar se- 



ries of simulations for a fixed value of p' = 29.93 (the same 
value as that for line 2)) and changing the value of Pcr 2 /e lh 
This p' is consistent with the thermodynamically stable point 
at fo-j/eo = 1.0 obtained above within the relative error of 
some percent. The given Pcr\jeo of lines 6) to 8) ranges 
Pct^/eq = 0.95 to 1.1. A snapshot of the system at 8) is 
presented in Fig. 2(b). Despite different intensive parameter 
sets, all the systems of 6) to 8) show the globally-anisotropic 
defect-free string-like assembly similar to Fig. 2(a). Long- 
lived defects are absent in these systems. All the lines 6) to 
8) in Fig. 3(a) are fluctuating around lines 1) to 5), i.e. the 
vicinity of per 2 =s 0.451. For the values of Pcr 2 /eo outside this 
range, N quickly diverges or vanishes. From these data, we 
recognize that the present thermodynamically stable point lies 
between 7): Pcr\lea = 0.98 and 8): Pa- 2 /e = 1.1. This ther- 
modynamically stable point is equal to the one we obtained 
above within a relative error of some percent, as is expected. 

Once the thermodynamically stable point is determined ac- 
curately in /iPr-ensemble, we can exchange the ensemble to a 
conventional one, e.g. A^fT-ensemble or A^V F-ensemble, and 
can perform a longer simulation run, which is free from the di- 
vergence or the vanishment of N. As is expected, we can per- 
form longer simulation runs even with the yuPF-ensemble if 
the thermodynamically stable point is determined with higher 
accuracy. 

Next, we check the stability of our simulation method. 
Starting from the instantaneous microstate shown in Fig. 2(a), 
we resume simulation runs after disconnecting the reservoirs 
1 or 2, which corresponds to simulations with NPT or pVT- 
ensemble. 

Snapshots of the system and the time evolution of per 2 in 
these resumed simulation runs are given in Figs. 4 and 5 re- 
spectively. Both of these simulation runs preserve the same 
string-like assembly. In A^Pr-ensemble case, a value of per 2 
that is close to that for three-reservoirs method is also ob- 
tained. Therefore, the equilibrium state obtained via three- 
reservoirs method does not change even after the ensemble is 
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FIG. 5. Time evolution of pa\ in simulations with NPT and fiVT- 
ensembles that are resumed from the instantaneous state shown in 
Fig. 2(a). In these simulation runs, the number of the particles is 
A' ~ 500 with minor fluctuations. For reference, the result of the 
simulation with /jfT-ensemble is also shown. 



the system size, in ^Vr-ensemble shown in Fig. 5. 

These simulation results show that the physical properties 
of the system, in conventional ensembles, sensitively depend 
on N and the system box size while they do not strongly de- 
pend on N in //Pr-ensemble. 



FIG. 4. Snapshots of the model polymer-grafted colloidal system at 
ksT/eo = 0.12. Black dots represent the centers of the particles and 
grey lines denote networks of overlaps between the particles. From 
the instantaneous state of Fig. 2(a), simulations are resumed after 
disconnecting the reservoir 1 or 2, i.e. simulations are resumed with 
A'ProryuVT-ensemble. (a): simulation result at 1.1 x 10 7 MCS, after 
resuming the simulation with A^T-ensemble. (b): simulation result 
at 1.6x 10 7 MCS, after resuming the simulation with /jVT-ensemble. 



switched. These results justify our method for determining 
the thermodynamically stable point with yt/Pr-ensemble. 

Here, one can recognize occurrence of a few defects and 
undulation of the string-like assemblies after switching the 
ensembles. This should be attributed to the difference in the 
nature of the thermal fluctuations of a finite system between 
different ensembles. For example, when a particle is removed 
from a perfect triangular crystal in //PF-ensemble and simul- 
taneously the ensemble is exchanged to A^PF-ensemble, long- 
lived defects are created in the crystal. Figure 4(a) shows this 
finite size effect. In the thermodynamic limit where the system 
size becomes infinitely large, such a difference should van- 
ish. We also recognize a slight drop in per 2 in //VF-ensemble 
by approximately 2%. This is also related to the finite size 
effect of the /iVF-ensemble, which will be further discussed 
in section HE 2. The undulation of the string-like assembly 
shown in Fig. 4(b) corresponds to zigzag instability typically 
observed in convection rolls in a fluid slab when its natural 
periodicity is suddenly changed 34 . This zigzag instability is 
consistent with the slight drop in pcr 2 v i.e. a slight increase in 



2. Alder transition of the outer cores. 

At extremely low temperature, the repulsive square-step of 
<p(r) becomes far higher than the thermal energy, kgT. Due to 
this extremely high potential energy barrier, the phase behav- 
ior of the system is almost identical to the behavior of hard 
particle systems with diameter cr 2 if the system volume V ex- 
ceeds the close-packed volume (area) of the outer cores of 
the particles, denoted by Vq (= No- 2 2 V3/2 in 2 dimensions) 2 . 
Therefore, crystalization of the hard particles at high density, 
called Alder transition 35 , occurs in our system at V/Vq m 1.3 
and low temperature 2 . Here we simulate this triangular crys- 
tal of the outer cores of our colloids with the diamter o"2 at 
low temperature ksT/eo = 0.1, where we fixed the parame- 
ters (fcr 2 /e = 0.45,//' = 18.42) that are determined via the 
iterative refinement of P and p.. The stability of this thermody- 
namically stable point is discussed later in the present section. 

Simulation with the yuPr-ensemble is started from an initial 
state with V/cr 2 = 4778.38 and N Q = 1089. We prepared this 
initial configuration by removing the particles from the equi- 
librium configuration with 33 x 38 = 1254 particles arranged 
on a homogeneous triangular lattice, which results in an inho- 
mogeneous particle configuration. 36 Such an inhomogeneous 
configuration is swiftly equilibrated in /vPF-ensemble as is 
shown in Fig. 6. A snapshot of the system simulated with the 
/ifT-ensemble, which is presented in Fig. 6, shows the defect- 
free triangular crystal. Temporary small defects due to ther- 
mal fluctuation are sometimes found in the system, whereas 
long-lived defects are absent. Time evolution of per 2 and 
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FIG. 8. Time evolution of pa\ in A'/T-ensemble at {k B T /eq 
0.1, Pa 2 ,/ eo = 0.45). 




FIG. 6. Snapshot of the model polymer-grafted colloidal system in 
fiPT -ensemble at {k B T/e = 0.1, Pa\/e = 0.45,//' = 18.42) and 
N = 1089 at 9 x 10 s MCS. The initial system volume is set at V/aj = 
4778.38. Black circles represent the inner cores of the particles and 
white ones denote the outer cores. 
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FIG. 7. Time evolution of pa\, (a), and N, (b) and (c), in the simula- 
tions with the /z/T-ensemble at (k B T/eo = 0.1, Pa 2 ,/eo = 0.45,//' = 
18.42). Black lines denote the results of that are started from the ini- 
tial condition with V/a\ = 4778.38 and N = 1089 and grey lines 
from that with V/a\ = 5647.18 and N = 1221. 



is plotted in Fig. 7. Small fluctuation in N indicates the high 
stability of this defect-free crystalline state. Different initial 
particle configuration also results in a similar defect-free crys- 
talline state with the same average per 2 . Time evolutions of 
per 2 and N for simulation on the system with different No and 
initial V/cr 2 are also plotted in Fig. 7, which show the same 
defect-free structure. This demonstrates the high stability of 
the defect-free crystalline state obtained via three-reservoirs 
method. 

Next, we compare physical characteristics of pPT- 
ensemble with these of the conventional ensembles. For this 



purpose, we perform simulations also with the conventional 
ensembles, i.e. NPT and //VF-ensembles, with the same pa- 
rameters (k B T/e = 0.1, Per 2 /e^ = 0.45,//' = 18.42). 

Different from the //fT-ensemble case, in the present NPT- 
ensemble case, we have to manually tune the value of N 
so that the perfect ordered equilibrium structure can be ob- 
tained. For this reason, we perform a series of simulations 
with A^Pr-ensemble for all the values of N within an interval 
1070 < N < 1089, where p v = l/N and 1 MCS = N sim- 
ulation steps are fixed and the initial system volume is still 
V/cr 2 = 4778.38. The only parameter that is changed from 
the above //Pr-ensemble simulation is N. Typical examples 
of the time evolution of per 2 in these simulation runs are plot- 
ted in Fig. 8. For reference, the result with the parameter op- 
timized in //Pr-ensemble, N = 1216, is also plotted in this 
figure. The system with this optimized parameter shows the 
defect-free triangular crystal and its value of per 2 is close to 
the results of the //fT-ensemble. However, with any value 
of N in the above interval, long-lived defects, mostly point 
defects, appear in the system, as is shown in Fig. 9. More- 
over, the average value of per 2 changes with the change in 
N slightly, nonmonotonically, and sensitively. This behav- 
ior is different from the results of the //fT-ensemble. This 
means that, although the external intensive variables (T, P) 
are fixed, physical properties of the system in A^fT-ensemble 
sensitively depend on the external extensive variable, N. 

In //Vr-ensemble, the system size has to manually be 
tuned. Actually, we try such a manual tuning by perform- 
ing the //Vr-ensemble simulation with various system sizes 
in an interval 3909.59 < V/cr\ < 5212.78. Aspect ratio of 
the system box is kept at a typical value for molecular simu- 
lation, L JL V = 1.00. Simulation parameters that are changed 
from those of the corresponding //Pr-ensemble simulation 
are V and N = 33 x 38 = 1254. We assume that 1 MCS 
= A^o simulation steps. All the systems we have simulated 
show the defect-free triangular crystals, since point defects 
are directly removed by the particle insertion and deletion 
processes. Long-lived defects are not found in these simu- 
lation runs with //Vr-ensemble. However, the average value 
of pcr\ sensitively and nonmonotonically depends on the value 
of V/crj. Typical examples of time evolution of per 2 are given 
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FIG. 9. Snapshots of the model polymer-grafted colloidal system simulated with ./VfT-ensemble at (k B T/eo = 0.1, Pct\Ieq = 0.45) and 2 x 10 7 
MCS. Black circles represent the inner cores of the particles and white ones denote the outer cores, (a): N = 1082. (b): N = 1089. 
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FIG. 10. Time evolution of pcr^ obtained with yuVT-ensemble at 
(k B T/€ = 0.1, /u' = 18.42). The black line denotes the result at 
V/o"j = 4900.94, which is the optimized value obtained with three- 
reservoirs method. 



abrupt and sensitive to the other parameters. However, in this 
/iVr-ensemble case, the point defects can rather easily be re- 
moved even in high-density states e.g. in a triangular crys- 
tal. On the other hand in iVPr-ensemble, p is changed by the 
change in the continuous dynamic variables, (L X ,L V ,L Z ). This 
results in a smaller change in the average p compared with that 
in juVr-ensemble. However, long-lived defects are frequently 
found in A^Pr-ensemble. Three-reservoirs method overcomes 
these disadvantages of fiVT and ATPF-ensembles and, at the 
same time, it inherits the advantages of these ensemble meth- 
ods. In /ifT-ensemble, both point defects and line defects can 
be removed easily, and the extensive variables can finely and 
spontaneously be tuned. 



F. Global equilibrium 



in Fig. 10. Although the intensive variables (T,/j) are spec- 
ified by the reservoirs, physical properties of the system in 
/iVr-ensemble significantly depend on the extensive variable 
V. The data for the simulation with Vjo\ = 4900.94, which is 
the optimized value obtained in the three-reservoirs method, 
are also plotted by a black line. We can confirm that the 
optimized value of per 2 is preserved during the simulation 
run. This result justifies the simulation results obtained using 
three-reservoirs method. 

In the Vr-ensemble simulations, the particle density p can 
be adjusted only discretely because of the discrete nature of 
N, and therefore the behavior of p in yuVr-ensemble is rather 



The simulation results in section II E show that it is a te- 
dious task to perform the manual tuning of N and the system 
size in conventional ensembles, whereas these extensive vari- 
ables are automatically and finely tuned in our yt/Pr-ensemble. 
In the conventional ensembles where at least one extensive 
variable is fixed, physical properties of the system are signif- 
icantly dependent on the value of such a fixed extensive vari- 
able. This illustrates that, for example, the equilibrium state 
of a system in A^fT-ensemble depends on N even though the 
intensive variables T and P are fixed. As another example, 
the equilibrium state of a system in A^Vr-ensemble changes 
with N and (L x , L y , L z ) even though the external intensive vari- 
ables T and p = N/V are fixed. The equilibrium state in con- 
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ventional ensembles is specified by the parameter sets of the 
external extensive variables as well as the external intensive 
variables. In the present article, we tentatively define the local 
equilibrium state as the equilibrium state at each parameter 
set of the external extensive variables with the fixed external 
intensive variables. 

In yuPr-ensemble, however, the extensive variables are 
finely and spontaneously tuned and the most stable state over 
the local equilibrium states at the given (T,P,/i) is automat- 
ically obtained. In the present article, we tentatively call 
this equilibrium state in yuPr-ensemble the global equilibrium 
state. 

As these dependences of the physical properties on the ex- 
ternal extensive variables in conventional ensembles are re- 
garded as a finite size effect, the three-reservoirs method is a 
technique to remove the finite size effect of the conventional 
ensemble methods. Any local equilibrium states are, in the 
thermodynamic limit, identical to the global equilibrium state. 



G. Entropy and free energy calculation in yu/T-ensemble 

In molecular simulations, physical quantities of the simula- 
tion system are defined through the ensemble average over the 
probability distribution of the microstates in the phase space. 
The evaluation of the free energy of such a system, however, 
is equivalent to the evaluation of its partition function. As 
the partition function is not an averaged quantity over the en- 
semble, its evaluation demands a high dimensional integral 
over the whole phase space, resulting in unrealistically large 
amounts of computational cost. 

Instead of evaluating the partition function directly, a 
derivative of the free energy with a control parameter, e.g. 
pressure in the canonical ensemble, is calculated for the sake 
of evaluating the free energy in standard molecular simula- 
tion 10 . This technique, which is usually called the thermo- 
dynamic integration, gives the free energy difference between 
two different thermodynamic states. 

However, with this technique one cannot go across a first- 
order phase transition. For example, when the two thermody- 
namic states are located in a solid phase and in a fluid phase 
respectively, a first-order transition occurs in the middle of the 
integration path. Due to possible hysteresis at this transition 
point, forward and backward integration paths between the 
two states in general gives different values for the free energy 
difference. This problem also affects the other free energy cal- 
culation techniques, e.g. histogram reweighting technique 37 ' 38 
and the method of expanded ensembles 11 . In order to over- 
come this difficulty, we need to find appropriate reference 
states, whose free energy values are already known 8 ' 10 , e.g. 
Einstein solid for the free energy calculation of crystals 6 ' 39 . 
These reference states provide reversible integration paths. If 
we cannot find such reference states, we have to find integra- 
tion paths that bypasses the first-order transition line, e.g. an 
integration path that is arranged with the help of artificially 
introduced external fields 5 ' 6 ' 40 ' 41 . After setting a reversible in- 
tegration path using these techniques, we run simulations at 
a large number of state points along the integration path, and 



integrate the derivative of the free energy along this path. In 
addition to the discretization error in the integration along the 
path, occurrence of defects in the system also affects the ac- 
curacy of the free energy evaluation of the ordered structures. 

Here we propose, based on Euler equation in thermody- 
namics, a convenient method for the free energy calculation 
using the systems connected to the three reservoirs. The en- 
tropy of the system per particle, denoted by s = S /N, satisfies 
Euler equation, 

s(T,P, f i) = l^-(PV- f iN+U + E k )\ , (15) 

where denotes the total kinetic energy of the system, U de- 
notes the total potential energy, and (• • • )t,p,h is the ensemble 
average at (T, P,/j). In the right-hand side of this equation, the 
3 intensive variables, T, P, and // of the system, relax to the 
equilibrium values that are equal to those given by the reser- 
voirs. {Ek/N)T,p,/j is determined via the equipartition theorem, 
e.g. (3/2)k B T for monatomic molecules. The ensemble av- 
erages {V/N)t,p,h and (U/N)t,pjj can directly and readily be 
evaluated through the simulation runs of the three-reservoirs 
method. Therefore, according to this Euler equation, s(T, P,[i) 
is determined from our simulation at one state point (T, P,/i), 
which means that this free energy evaluation method requires 
far smaller amounts of computation than the other standard 
methods do. The free energy of the system, e.g. Helmholtz 
free energy, Gibbs free energy, and grand potential, is also ob- 
tained from this result in a similar manner. Since our entropy 
and free energy calculation method is free from any thermo- 
dynamic integration paths, the first-order phase transition does 
not affect our evaluation method. In addition, with the three- 
reservoirs method we can easily eliminate the defects which 
is the main origin of the error in the free energy evaluation for 
the ordered phases. This raises the accuracy of the evaluation. 

In our molecular simulation, fi includes the thermal de 
Broglie wave length A, as is discussed in appendix A, i.e. 
eq. (Al). h and mm A cancel when the free energy difference 
between two different state points is calculated. Therefore, we 
do not have to take care of these h and m. 

Although the above method based on Euler equation is, in 
principle, applicable to simulations with the other ensembles, 
e.g. NVT and AffT-ensembles, one has to measure intensive 
variables P and/or /j. as ensemble averaged values. Such a 
procedure requires a computationally expensive analysis. For 
example, the measurement of fi, i.e. Gibbs free energy per 
particle, demands a vast amount of simulation, especially in 
high particle density regions 5 ' 6 ' 42-46 . However, in the \xVT- 
ensemble or /iPF-ensemble, we do not have to evaluate /j. be- 
cause it is already given by the reservoir. Same is true for the 
evaluation of the pressure P. As a result, with the use of the 
/iPr-ensemble, we can skip tedious evaluations of the inten- 
sive variables because all the essential intensive variables yu, 
P, and T are already specified by the reservoirs. 

In addition, metastable structures and defects, which fre- 
quently appear in conventional ensembles, affect the results of 
this entropy and free energy evaluation. When the system is 
in metastable states or outside the global equilibrium state, the 
intensive variables given from the reservoirs are inconsistent 
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with the values of these variables in the simulation system. 
This means that both fi and P need to be analyzed in the sim- 
ulation rather than to use the specified value by the reservoirs. 
However, as the defects are quickly eliminated and the sys- 
tem reaches the global equilibrium in the simulations with the 
/iPr-ensemble, the evaluation of the free energy is free from 
the above problem associated with the metastable states and 
the local equilibrium states. 

1. Construction of equilibrium phase diagrams. 

When one tries to construct the equilibrium phase diagrams 
using conventional simulation methods, candidates for the 
equilibrium structure have to be chosen prior to the simula- 
tion and the free energy of each candidate should be measured 
and compared with each other with high accuracy, for exam- 
ple, with a typical error level 5 ' 6 ' 28 of 10~ 4 to 10~ 6 . A variety 
of phases, e.g. crystals, the string-like assemblies, and the 
other ordered and disordered structures should be considered 
the candidates for the equilibrium phase. In actual calculation, 
we empirically select some of these potential candidates and 
discard the others. However, we cannot deny the possibility 
that we have discarded the true equilibrium structure in this 
selection process. In addition, there is another possibility that 
the equilibrium structure is a totally new structure which has 
not been discovered yet. 47 

With the use of //PF-ensemble, however, one can obtain 
the global equilibrium structure directly as a result of the 
fine tuning of the extensive variables. Therefore, the above- 
mentioned problem encountered in the construction of the 
phase diagram using the standard ensembles can be avoided 
when we use /iPr-ensemble. 

III. CONCLUSIONS 

We have studied thermodynamics, statistical mechanics, 
and molecular MC simulation algorithms of yuPr-ensemble. 
Guggenheim formally introduced Boltzmann factor (statisti- 
cal weight) of this ensemble with an assumption that the av- 
erages of extensive variables can be determined 13 . However, 
this assumption contradicts the indetermination of extensive 
variables in //PF-ensemble. In addition, other characteristics 
of this ensemble were totally absent in Guggenheim's dis- 
cussion and other early works 14-16 . These early researchers 
concentrated on the formalism, i.e. mathematical aspects of 
the partition function of this ensemble, only at the thermo- 
dynamically stable point 17 . Physical aspects of the ensemble 
have not seriously been discussed. In the present work, we 
have shed light on these problems and have discovered ther- 
modynamic and statistical mechanical characteristics of iiPT- 
ensemble, e.g. the thermodynamically stable point, thermody- 
namic degrees of freedom, thermodynamics outside the ther- 
modynamically stable point, maximization of entropy density, 
quick equilibration due to shortcuts in the phase space, etc. 
We have also shown that //PF-ensemble is built as a com- 
bination of the 3 ensembles each of which is combined to 2 



reservoirs, i.e. NPT, p.VT, and p.PS -ensembles. The 3 cor- 
responding free energy densities, i.e. Gibbs free energy per 
particle, grand potential per volume, and (E - fiN + PV)/S, 
are simultaneously minimized in yuPr-ensemble, rather than 
the thermodynamic potential of yuPr-ensemble, denoted by 
O. 

We have proposed a molecular MC simulation method 
based on /iPr-ensemble (three-reservoirs method) 12 . We can 
show that this three-reservoirs method gives a physically ac- 
ceptable ensemble, which allows us to trace the physical tra- 
jectories in the phase space. Since three-reservoirs method 
is built as a combination of conventional NPT and fiVT- 
ensembles, programming is lighter than other advanced tech- 
niques. In addition, the thermodynamically stable point is de- 
termined according to Gibbs-Duhem equation in a short com- 
putational time. These features mean that three-reservoirs 
method requires a small amount of preparation and that we 
can quickly start production simulation runs, although other 
advanced techniques demand a large quantity of complicated 
preparation, e.g. advanced programming, the precise adjust- 
ments of the artificial weights necessary for multicanonical 
technique, and accurate free energy measurement essential for 
the expanded ensemble technique 11 . 

These advantages over other simulation techniques facil- 
itate and reduce the total work flow of our three-reservoirs 
method compared with the conventional methods such as 
NPT, fiVT, or multicanonical method. Furthermore, only 
with the three-reservoirs method, we can simultaneously and 
automatically tune the number of particles N and the system 
size to obtain the equilibrium ordered state. This unique ad- 
vantage of the three-reservoirs method could enhance the un- 
derstanding of those systems that were obtained via the other 
standard simulation techniques. For example, for perfect crys- 
tals, both N and (L x ,L y ,L^j must be integer multiples of the 
unit structure of the crystal structure, which is in general not 
known a priori. In principle, by measuring the free energy 
density of various structures at each set of N and (l x , L y , L z ), 
these extensive parameters can manually be tuned in simu- 
lation of ATPr-ensemble, //VF-ensemble, multicanonical en- 
semble, etc. In practice, however, this manual tuning requires 
much computational effort. On the other hand, with the use 
of our three-reservoirs method, we can automatically achieve 
such optimization. 

For a solid at finite temperature T, the system with the fiPT- 
ensemble reaches a globally-anisotropic defect-free ordered 
state as the equilibrium state by crossing the metastable states 
through the shortcuts in the phase space due to the additional 
degrees of freedom. On the other hand, in conventional en- 
sembles, physical properties of the system sensitively and dis- 
cretely depend on N and/or (l x ,L x ,L^ even though external 
intensive variables are fixed. This results in a requirement of 
the manual tuning of these extensive variables to obtain the 
global equilibrium state. 

These results illustrate that our method can be applied to a 
variety of physical systems for the sake of studying ordered 
structures in equilibrium at finite T, e.g. lamellae composed 
of diblock copolymers, smectic phase of liquid crystals, and 
fluid bilayer membranes. Three-reservoirs method can also 



13 



be applied to numerical calculation of the equation of state 
that relates p, P, and T. We have also shown that the entropy 
and the free energy can quickly be evaluated in pPT -ensemble 
based on Euler equation. This feature is essentially important 
in the construction of the phase diagram of condensed materi- 
als. 
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In simulation with the /iVr-ensemble 8 ' 10 ' 50 , particles are 
inserted into and deleted from the system in addition to 
Metropolis trial displacement of particles. In one simulation 
for this ensemble, these steps are included in: 

i) with probability pc/2, trial particle insertion into the sys- 
tem 

ii) with probability pc/2, trial particle deletion from the sys- 
tem 

iii) with probability 1 - p G , trial displacement based on 
Metropolis algorithm, i.e. perturbation to one particle 

is chosen, where p c is a constant fixed in an interval < pc < 
1 . Algorithms of trial particle insertion and deletion are dis- 
cussed in the following. 



Appendix A: Molecular MC technique in fiVT and 
A'/T-ensembles 



In this appendix, we show molecular MC simulation tech- 
nique 8 ' 10 ' 48 ' 49 for single component systems based on the pVT 
and Wr-ensembles. We use the following notations: r, de- 
notes the spatial coordinates of the particle ;', m mass of a par- 
ticle, kgT the thermal energy, U the potential energy of the 
system, and h Planck's constant. The chemical potential of a 
system consisting of real particles in the canonical ensemble 
(iVVr-ensemble) is given by, 



Pka(T,p) 



Pidenl(T,p) /£exccs 



„(r,p) 



where 



and 



/"ideal (T>P) 



knT 



p:=N/V, (Al) 



k B T 



:=log(/l 3 p), 



(A2) 



a. Particle insertion. 

We assume that one particle is inserted into the system that 
is composed of N particles. The position of this inserted par- 
ticle r;v + i is chosen uniformly over the system box. The co- 
ordinates of the N particles, (r\ , . . . , are fixed during this 
particle insertion step. The state after the insertion, i.e. the 
state of N + 1 particles, is accepted as a new state with the 
probability, 



acc(JV -» N + 1) 
= min 1 , 



V 



A 3 (N + 1) 



exp 



£ Y excess 



k B T 



(A4) 



U™ cess :=U(r u ...,r N+1 )-U(ru...,r N ). 



A function min(ai,fl2) returns the smaller of two arguments, 
«i and «2. If the trial insertion is rejected, the state before the 
insertion is kept for the next simulation step. 



A := 



\l2nmk B T 



(A3) 



Here, Pife a \(T,p) denotes the chemical potential of an ideal gas 
in A^Vr-ensemble and A is called the thermal de Broglie wave 
length. The quantity p excess (T,p) denotes the excess chemical 
potential that is originated from the interaction between the 
real particles 6 ' 10 . 

Simulation results in /iVr-ensemble are independent of A. 
This means that A only appears in the chemical potential p 
at the reference point in /iVr-ensemble simulation, which is 
discussed in section A 1. 



1 . MC simulation method in fi VT-ensemble 

In the present section, we discuss the MC method of a sin- 
gle component system in the grand canonical ensemble, which 
is also called yuVr-ensemble because p, V, and T are kept 
fixed. 



b. Particle deletion. 

We assume that one particle is randomly chosen and at- 
tempts to be removed from the present system composed of 
N particles. This chosen particle, denoted by index is re- 
moved from the system with the probability, 



acc(N -> N - I) = min 



AN 



1 



' V excf l -U AA + -2-1 

V CA P[ /;g7" excess kgT J / 



(A5) 



taxless := U (n, . . . ,r N ) - U (n, . . . , rj-u r j+u . ..,r N ). 

If the trial deletion is rejected, the state before the deletion is 
kept for the next simulation step. 

When p le . d \ is substituted for p in eqs. (A4) and (A5), these 
acceptance criteria are independent of A. This illustrates that 
simulation results in //VF-ensemble are free from the actual 
value of A. Therefore, A appears only in the chemical poten- 
tial p at the reference point in the yuVr-ensemble simulation 6 . 
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2. MC simulation method in A7 J 7 -ensemble 

MC simulation method in isobaric-isothermal ensem- 
ble 10 ' 51 ~ 53 , also called A/PF-ensemble, is briefly summarized 
in this section. In addition to N and T, the system pressure, 
denoted by P, is given from the outside, in this ensemble. 

In simulation of A/PF-ensemble, the system size is 
changed, in addition to the trial move of particles. In one sim- 
ulation step, we include the following steps: 

i) with probability p v , trial system size change 

ii) with probability 1 - p v , trial displacement based on 
Metropolis algorithm, i.e. perturbation of one particle 

is chosen, where p v is a constant fixed in an interval < p v < 
1 . Algorithms of trial system size change are discussed in the 
following. 



a. Trial system size change in NPT-ensemble. 

To change the system box size, we use the following 4-step 
algorithm. 

We assume that the system box size, denoted by (l x , L y , L z ), 

is changed to a new system size, (l' x , L' y ,L' z J. Unlike the 
conventional MC simulations in A/TT-ensemble based on 
McDonald's method 10 , each element of the system size 
(l x , L y , L z ) is independently changed in our algorithm. 

1) The new system size is chosen, 



L' X = L X +AL(\-2^ X ), 
L y =Ly+AL(l-2t y ), 
L' Z =L Z +AL(\-2^ Z ), 



(A6) 



where AL is a small constant length, g x ,jj y , and £. are 
random numbers uniformly distributed over an interval 
< ^ 1. and V = L' x L' y L' z denotes the volume 

of the new system box. 

2) Coordinates of all the particles before the system size 
change, r, = (x,-, y,-, z ( ), are homogeneously scaled, 

r\ = {{L' x IL x ) Xi , {L' y /L y ) yi , (L' z /L z ) Zl ) , 1 < i < N. (A7) 
This is new particle coordinates after the change. 

3) The potential energies of the systems before and after the 
system size change, U(r\,..., r N ) and u(r\,..., r' N ^j re- 
spectively, are calculated. 

4) The new system size and particle coordinates are accepted 
with the probability, 



acc(V -> V) = 
min 1 , exp 



k R T 



U(r' v ...,r' N )-U(n,...,r N ) 



+P(V - V) -Nk B T log 



V 



(A8) 



If this trial system size change is rejected, the state before 
the trial is kept for the next simulation step. 



Appendix B: Ensemble average at each A' and V in 
yuPr-ensemble 

The ensemble average at each N in //Pr-ensemble is equiv- 
alent to the average in A/PF-ensemble at the same N. We il- 
lustrate this in the present section. 

We assume that a finitely long simulation run is performed 
in /iPr-ensemble at a thermodynamically stable point. Af- 
ter the equilibration of the simulation system, a finitely large 
number of microstates of the system are visited. Unnormal- 
ized Boltzmann factor of the system is denoted by w,, where 
the suffix i represents these microstates. Si m is Kronecker 
delta. /?= l/k B T. 

After the equilibration of the system, the ensemble average 
of an intensive physical quantity, A, obtained in the present 
simulation run is: 



(A) 



N=0 XiWi ^ X,<W W « 
hi 2* S Ni ,N ^fm exp [-p {PVi - l iN+ Ud] 



N=0 



ZiS N „NAi exp [-13 (PVi + Ud] 



Zi6 NuN ex P [-/3(PVi + Ud] 
= Y,f(N)(A) T ,p,N, 



(Bl) 



N=0 



where the summation 2/ mns ove r the microstates visited in 
the simulation, and f(N) is defined as, 



/(AO := 



(B2) 



and (A)tj>,n denotes the ensemble average of A in NPT- 
ensemble. f(N) denotes the occurrence probability of N in the 
simulation run. An expression similar to eq. (Bl) also holds 
for the ensemble average of A in /iVr-ensemble. This result 
illustrates that the short-time average of A in //PF-ensemble is 
approximated by the ensemble average in A/PF-ensemble or 
/iVr-ensemble. 

When the finite size effect of the system, which has been 
discussed in sections II E and II F, is small, {A) T ^ N is depen- 
dent on (T, P) and independent of N. Therefore, eq. (Bl) 
yields {A) TPll = {A) TPN = A(T,P). A similar relation, 
(A)t,Pjj = (A)t,h,v = A(T,n), also holds, where (A) T ^ V de- 
notes the ensemble average of A in //VF-ensemble. Therefore, 
(A) T ^=A(T,P)=A(T,n). 
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